set seed 23423

local colist HU SE RO PT PL NO IT GB FR ES DK DE CZ BG AT

** loop over the NUTS definitions and productivity variables
foreach nuts in 2 3{
foreach ltfpvar in w1_lvapw  w1_log_TFPr_tfas{
foreach co in  `colist'{
	di "Country is `co'"
	quietly use "$BvD_input/Amadeus_Europe_2010_2019_coded_latest_NUTS`nuts'_`co'.dta", clear
	quietly gen order=_n
	
	*Shorthand for the productivity variables we loop over: "lA"
	quietly rename `ltfpvar' lA
	quietly gen N_total=_N	
		
	*Generate group IDs
	quietly egen JT=group(J T)
	quietly egen ST=group(S T)
	quietly egen SJT=group(S J T)
	quietly egen tag_SJT=tag(SJT)
	quietly egen tag_ST=tag(ST)
	
	*Counts (of firms and cells)	
	quietly egen count_xi=total(tag_ST)
	quietly egen count_tau=total(tag_SJT)
	quietly gen count_firms=N_total
	quietly gen countpercell_xi=count_firms/count_xi
	quietly gen countpercell_tau=count_firms/count_tau
	
	*Shares of industry-location cell in location cell 
	quietly bys SJT: egen emp_SJT=total(empl) if lA!=.
	quietly bys ST: egen emp_ST=total(empl) if lA!=.
	quietly gen empshare_SJT=emp_SJT/emp_ST
	
	
	*Generate location averages of productivity: FEs from firm-level data

		*Location-industry cell average ("tau")
			* Step 1: generate industry-location FEs, weighted by firm emp
		quietly reghdfe lA [aw=empl], absorb(temp_tau=SJT) keepsingletons
			* Step 2: demean those FEs (unweighted cells regression; recall that tag_SJT is =1 just for one firm per cell)
		quietly reghdfe temp_tau if tag_SJT, absorb(JT) resid keepsingletons
			* Step 3: the residual from step 2 is the "tau" in the paper (only assign to tag-firms)
		quietly predict tau if e(sample), resid	
			
		*Location cell average ("xi")
			* Step 1: take taus from above but now average at location level with 
			* employment shares at industry-location level (recall that tau is 
			* just for 1 firm per cell, so this is a cell-level regression)
		quietly reghdfe tau [aw=empshare_SJT], absorb(xi_temp=ST)  keepsingletons
			* Step 2: center around 0 (recall we only have one year, set T=1)
		quietly reghdfe xi_temp if tag_ST, absorb(T) resid keepsingletons
			* Step 3: the residual from step 2 is the "xi" in the paper (only assign to tag-firms)
		quietly predict xi if e(sample), resid

		*Calculate variances of xi, save var and N
		quietly tabstat xi if tag_ST, stat(N var mean) save
		matrix stats		= r(StatTotal)
		quietly gen N_xi	= stats[1,1]
		quietly gen var_xi= stats[2,1]
		
		*Calculate variances of tau (see paper Section 3 for discussion of weighting), save var and N
		quietly tabstat  tau [aw=empshare_SJT]  if tag_SJT, stat(N var mean) save
		matrix stats		= r(StatTotal)
		quietly gen N_tau	= stats[1,1]
		quietly gen var_tau	= stats[2,1]	
			
			
		*Cacluate p90-p10 log differences
		foreach var in tau xi {
		*generate weights for tau and xi
			tempvar weight
			if("`svar'"=="tau"){
				gen `weight'=empshare_SJT
				local tag SJT
				}
			else{
				gen `weight'=1
				local tag ST
				}
		_pctile `var' [aw=`weight'] if tag_`tag',  p(10, 90)
		local `var'_p10 =r(r1)
		local `var'_p90=r(r2)
		quietly gen `var'_p9010=``var'_p90'-``var'_p10'
		}
		
		*keep country-level information only (keep first line of code)
		quietly keep if _n==1
		quietly keep  var* count* *_p9010
		quietly gen co="`co'"
		if("`co'"=="HU"){
			quietly save "$BvD_input/BvD_results_`ltfpvar'_`nuts'.dta", replace
			}
		else{
			quietly append using "$BvD_input/BvD_results_`ltfpvar'_`nuts'.dta"
			quietly save "$BvD_input/BvD_results_`ltfpvar'_`nuts'.dta", replace
			}		
	}
quietly use "$BvD_input/BvD_results_`ltfpvar'_`nuts'.dta", replace

* generate final data set at country-level (separate files for NUTSxProductivity)
gen tfpvar="`ltfpvar'"
if ("`ltfpvar'"=="w1_lvapw"){
	quietly save "$BvD_output/BvD_results_`nuts'.dta", replace
	}
	else{
	quietly append using "$BvD_output/BvD_results_`nuts'.dta"
	quietly save "$BvD_output/BvD_results_`nuts'.dta", replace
	}
 }
}
